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The distribution of counterions and the electrostatic interaction between two similarly charged 
dielectric slabs is studied in the strong coupling limit. Dielectric inhomogeneities and discreteness 
of charge on the slabs have been taken into account. It is found that the amount of dielectric 
constant difference between the slabs and the environment, and the discreteness of charge on the 
slabs have opposing effects on the equilibrium distribution of the counterions. At small inter-slab 
separations, increasing the amount of dielectric constant difference increases the tendency of the 
counterions toward the middle of the intersurface space between the slabs and the discreteness of 
charge pushes them to the surfaces of the slabs. In the limit of point charges, independent of 
the strength of dielectric inhomogeneity, counterions distribute near the surfaces of the slabs. The 
interaction between the slabs is attractive at low temperatures and its strength increases with the 
dielectric constant difference. At room temperature, the slabs may completely attract each other, 
reach to an equilibrium separation or have two equilibrium separations with a barrier in between, 
depending on the system parameters. 



PACS numbers: 87.10.-e, 87.15.-v, 82.70.-y 
INTRODUCTION 

Electrostatic interactions are of great importance in 
many biological and soft matter systems. For example, 
electrostatics plays a key role in proteins structure pQ, 
interaction of DNA and proteins with charged ligands 
[2], DNA packaging and condensation [3 , conformation 
of polyelectrolytes in solution and their aggregation be- 
havior [H|S] and interaction between charged colloids and 
their collective behavior in solution [BJ . Numerous exper- 
imental and theoretical studies have been performed to 
understand these phenomena and to elucidate the role of 
electrostatic interactions. The study of the counterions 
distribution and electric potential around DNA [7J, the 
interaction between two DNAs [5j and the same studies 
on similarly or oppositely charged surfaces [5J are exam- 
ples of these studies. 

The difference between dielectric constant of water 
(~ 80) and that of organic materials (~ 2), silica (~ 4) 
and air (~ 1) is an important factor in electrostatic in- 
teractions of related systems. The role of such dielectric 
inhomogeneity has been investigated in systems like ionic 
channels [TO], ions or colloids at the air- water interface 
[TT] , copolymers in electric field [T^] and polyelectrolytes 
adsorbed to an interface [T3]- More efforts are needed 
to investigate the effects arisen from dielectric inhomo- 
geneity on distribution of the counterions around charged 
objects and the interactions between them [l4l H"5] . 

When the electrostatic interactions in a system are 
much stronger than the thermal energy, mean field the- 
ories such as Poisson-Boltzmann formalism are not suit- 
able for describing the behavior of the system. In this 
situation, known as the strong coupling (SC) limit [IB] . 



another theory, which is derived from the virial ex- 
pansion of the system partition function, can be used. 
This theory has been used to study electrostatic interac- 
tions between charged objects and counterion distribu- 
tion around them. Such studies have been performed 
for two parallel charged plates considering their sur- 
face charge density as uniform [17] and discrete [HI [19] . 
Also, electrostatic interaction between two overally neu- 
tral plates of randomly distributed annealed or quenched 
charges of both signs has been studied [50J . Corrections 
to the theory of Ref. [TBJ are suggested to describe the 
system behavior in a broader range of the parameters, 
namely the intermediate coupling regime [21 . Such cor- 
rections are recently formalized by expansion of the par- 
tition function around the ground state of the system in 
which the Wigner crystal is formed by the counterions 
[22] . This approach gives a correction to the prefactor of 
the second term of the theory of Ref. [16] to agree with 
Monte Carlo simulations. 

The system of two uniformly charged dielectric slabs 
has been studied taking into consideration the effects 
arising from dielectric inhomogeneities [HI H5] . It has 
been shown that in the SC limit the dielectric inhomo- 
geneities cause accumulation of the counterions in the 
middle of the intersurface space between the slabs. Also, 
it has been shown in these studies that the dielectric inho- 
mogeneities increase the repulsive pressure between the 
slabs at small inter-slab separations. 

A widely used simplification in the studies of charged 
systems is the consideration of the charge distribution on 
the surfaces of charged objects as uniform. Clearly, this 
approximation is feasible only when the ions in the solu- 
tion are far from the charged surfaces and also the sepa- 
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ration between charged surfaces are quite larger than the 
typical separation between charged residues on the sur- 
faces. When a point charge is very close to a big charged 
object or when two charged objects are very close to each 
other, discreteness of charge should be taken into account 
in calculation of the electrostatic interactions. In the SC 
limit, it has been shown that the distribution of countc- 
rions in the vicinity of a discretely charged plate is very 
different from that of a uniformly charged plate. In the 
former case, counterions have a strong lateral correla- 
tion with the surface charges and counterion density on 
the surface is much higher relative to the case of uni- 
formly charged surface [18] . For the case of two discretely 
charged surfaces of the same sign, most of the counterions 
crowd near the surfaces and the attraction between the 
surfaces is stronger compared to the case of uniformly 
charged surfaces [T5]. The same differences have also 
been observed in the weak coupling limit [23l [24] . 

The purpose of this paper is to consider the effects of 
both dielectric inhomogeneity and discreteness of charge 
on the electrostatic interaction between two similarly 
charged parallel dielectric slabs. We consider the charge 
distribution on the surface of each slab as square shaped 
uniformly charged regions of side length I whose centers 
are on a square lattice of spacing a, as is shown in Fig. [T] 
The case of I = a corresponds to uniform charge distribu- 
tion on the slabs and I — > corresponds to distribution 
of point charges. We use the Green function method to 
calculate interaction energies of each counterion with its 
images and the surface charges, and use it to obtain the 
counterion density profile in the SC limit. Despite pre- 
vious works (e.g. Refs. [IH112I]), all the Fourier modes 
of the charge distribution are taken into consideration 
here. This allows us to study the extreme limit of point 
charges on the surfaces, where the system behavior is 
noticeably different (see below) . Interaction between the 
slabs is calculated by two approaches: The first approach 
is valid for low temperatures and the point-charge distri- 
bution on the slabs in which interactions of the counte- 
rions with each other is taken into account. The other 
approach is valid for room temperature and arbitrary val- 
ues of I, which is accurate in the SC limit where the effect 
of counterions interactions on the pressure between the 
slabs can be ignored. 

We find that consideration of the dielectric inhomo- 
geneity in electrostatic interactions of the system in- 
creases the tendency of the counterions to gather in the 
middle of the intersurface space between the two slabs. 
Taking into account the discreteness of charge on the 
slabs however, increases the counterion density in the 
vicinity of the slabs surfaces. Resulted from these two 
competing effects, when the slabs are in far or close sep- 
arations the counterions are mostly distributed in the 
vicinity of the surfaces or in the middle of intersurface 
space between them, respectively. In the limit of point- 
charge distribution of the surfaces, independent of the 



dielectric inhomogeneity strength, density of the coun- 
terions is maximum in the vicinity of the slabs surfaces 
at all inter-slab separations. It is found that the inter- 
action between the slabs, at low temperatures is always 
attractive and its strength increases with increasing the 
difference between dielectric constants of the slabs and 
that of the environment. At room temperature, depend- 
ing on the strength of charge discreteness and dielectric 
inhomogeneity, the system is found to have two differ- 
ent behaviors: the slabs stand in an equilibrium distance 
from each other or completely attract each other after 
passing a repulsive barrier in the intermediate separa- 
tions. It is also found that for a range of the system pa- 
rameters there are two equilibrium separations between 
the slabs corresponding to two minima in the free energy 
of the system. In the absence of dielectric inhomogeneity, 
by changing the charge distribution on the surfaces from 
uniform to point charges, the entropic repulsion barrier 
between the surfaces disappears and they attract each 
other at all values of inter-slab separation. These results 
are not observed in previous works and show the im- 
portance of considering the discreteness of charge. The 
effects of the dielectric inhomogeneity and the charge dis- 
creteness on the pressure between the slabs are also stud- 
ied. 

THE MODEL AND THE GREEN FUNCTION 

We consider a system of two infinite dielectric slabs 
of thicknesses b and c with dielectric constant £2. The 
slabs are parallel and placed in a medium of dielectric 
constant e\ at separation 2D from each other (Fig. [I]). 
We assume that the charge on the inner surface of the 
slabs is distributed as a square lattice of square shaped 
regions of side / and charge density a. The centers of the 
square shaped regions form a square lattice of spacing a, 
as shown in Fig. [l] We define the cartesian coordinate 
system, in such a way that the z-axis is perpendicular to 
the surfaces and the origin is located in the middle of the 
separation area along the z direction. By this definition, 
the position of the charged surfaces of the slabs along the 
z direction is z* = ±D. The charge density of the slabs 
in three dimensions can be written as 

p(x, y, z) = aS(z — z*) x 

+00 

Q(x + ma + l/2)@(-x-ma + l/2) x 

771 — — OO 
+ OO 

®{y + na + l/2)e(-y-na + l/2), (1) 

n— — 00 

where S(z) is the Delta function and Q(x) is the theta 
function that is zero for x < and one otherwise. Keep- 
ing the value of al 2 fixed, I — > corresponds to the dis- 
tribution of point charges on a square lattice of spacing 



3 



a on the slabs surfaces, whereas I — a corresponds to 
the uniformly charged surfaces of the slabs. The case of 
< I < a describes the distribution of uniformly charged 
domains on the slabs surfaces. It is worth mentioning 
that to avoid any undesired singularities, a cut-off layer 
of thickness t is supposed on each slab where the coun- 
terions cannot penetrate. 
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We consider the system in the SC limit. The cou- 
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pling parameter is defined as S = fj- — in which / b = 

2 

47nTs~ir~T * s t ne Bjerrum length outside the slabs, and 

2 

jU = 2-nqi B< j s 1S ^ e Gouy- Chapman length, where as de- 
notes the surface charge density and is given by as = 
(al 2 )/a 2 . We note that q is the counterions charge and 
a is the surface charge density of the charged regions. 



b 2D c 

FIG. 1: (Right panel) Two dielectric slabs of dielectric con- 
stant £2 placed in a medium of dielectric constant ei. The 
slabs are of thicknesses b and c and their distance is 2D. The 
counterions are not permitted to enter a region of thickness t 
on the slabs. (Left panel) On the inner surface of the slabs, 
charges are distributed on a lattice of constant a. On each lat- 
tice site, there is a square charged region of side I and surface 
density a. 



To calculate the electrostatic potential and interaction 
energies of the system, we use the Green function method 
and solve the equation V 2 Gi(r, r') = — AnS (r — r') <5 i3 
with appropriate boundary conditions, G; = Gj+i and 
where i — I, 2 . . . , 5 represents the 
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regions shown in Fig. [T] and Si is the dielectric constant 
of region i. By defining A = §^pf|j the Green function 
in the region between the slabs is found as 



G 3 (r,r ) 



f 

2^ 



+ 00 



(e~ kz > + W(kc)e- 2kD+kz >){e kz < + W{kb) e - 2kD - kz < 
k(l - W(kc)W(kb)e' 4kD ) 



+ika;(x-x')+iky(y-y') 



(2) 



r 



where z > (z<) is the larger (smaller) of z and z', and the 
function W(x) is defined as 



-2x 



W(x) = A 



1 - A 2 e- 2x 



(3) 



To find the Green function in the presence of a single 
dielectric slab, one should substitute c = in Eq. [5] Now 
the electrostatic energy of the system can be calculated 
numerically using the calculated Green function and all 
the Fourier modes of the charge distributions. 



Distribution of the counterions 

In the SC limit, distribution of the counterions can be 
found by calculation of the electrostatic energy of the 
system. Contribution of a countcrion to the electrostatic 
energy of the system can be written as E = UA,int+U s i a b, 



where UA,int is the image-counterion interaction arising 
from the dielectric discontinuity across the interface, and 
U s iab accounts for the counterion interaction with the 
physical charges on the slabs. In order to find the image- 
counterion interaction energy of each counterion, we as- 
sume a point charge q at position fo = (xo, yo, zq) located 
between the two slabs. One can find the potential at po- 
sition r as 0a(f) = j^G3(r,ro), where the Green func- 
tion is introduced in Eq. Substituting this potential 
in U = h J p (r) <fi (r) dV, the total electrostatic energy 
of the system in the presence of one counterion can be 
obtained. This energy also contains the self-energy of 
the point charge, which is infinite and independent of 
the slabs parameters. The self-energy can be obtained at 
the limit of D — >• oo. After subtracting the self-energy of 
the point charge, the total image-counterion electrostatic 
energy of a counterion between the two dielectric slabs is 
found as 
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+ °° „ W(kc)e 2kz ° + W(kb)er 2kz ° + 2W(kc)W(kb)er 2kD 
dk 



2 Jo 



^+2kD 



W(kc)W(kb)e- 



2kD 



(4) 



The electrostatic interaction energy of a counterion trie slab containing square-shaped charged regions on its 
with the charges of the slabs, U s i a b, can be calculated surface can be found using the Green function of Eq. ^ 
using the relation U — q<p. The potential of one dielec- by substituting c = as 



?lslab 



E 



8a 
e 2 a 2 

2al 

s 2 a 2 

al 2 

2e 2 a 2 ^ 



+°° a -kz 



m,n=l 
+°° „-k 



(W(kb) + 1) 



x sin(k x l/2) sm(k y l/2) cos(k x x) cos(k y y) 



J2 Z(W ^ h) + l) x sin(ifc'/ L>)(u»s( //.,-) com////)! 



(5) 



where k x = 2™, k v = 2™, k = ^Vn 2 + m 2 and k' = 
2j ^ k . We denote the relative displacement of the two 
dielectric slabs containing square-shaped charged regions 
in the x and y directions by 4> x and <p y , respectively. 
Clearly, in addition to above mentioned parameters of 



the slabs, the interaction energy of the system depends 
on the values of cf> x and <f> y as well. Using the Green 
function of Eq. ([2]), the potential of the two dielectric 
slabs in the space between them can be written as 



- kz+kD + W(kb)e- kD+kz )(l + W(kb)) . „ , . „ , ., , „ , 
kk x k v (e 2kD ~W 2 (kb)e- 2k D) sm (M/2) sm(M/ 2 ) cos(k x x) cos(k y y) + 
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^M sinfe | /2)sin( y /2) 
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The ion interaction energy with the charges on the slabs 
is obtained using this potential and the relation U s i a b — 



,J 2slab- 



In the SC limit, S> 1, the counterions distribution in 
the system can be written as a virial expansion, and to 
the leading order one has |16| 



PSC = ae 



(7) 



where E = f/A.ini + U s i a b is the total electrostatic en- 
ergy of each counterion discussed above and a is a nor- 
malization prefactor that can be calculated from relation 
J pdV — N, which yields the total number of the counte- 
rions in the system. In the following sections, we use this 
relation and the calculated energy of each counterion to 
find the counterion distribution in the presence of one or 
two dielectric slabs of discrete surface charge. 
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FIG. 2: (color online) The total electrostatic energy per coun- 
tcrion for (a) A = 0.95 and (b) A = 0.1. In both plots, the 
solid blue lines correspond to l/a = 0.8, the red dashed lines 
correspond to l/a = 0.5, and the black dashed-dotted lines 
correspond to l/a = 0.2. 



One dielectric slab with discrete surface charge density 

As mentioned above, the electrostatic energy of a coun- 
terion of charge q in the presence of one dielectric slab 
can be written as 



E 



nslab 



+ 00 



dkW(kb)e 



-2k z 



(8) 



where the first term is the electrostatic energy of the 
counterion with the physical charges on the slab and the 
second term shows the image-counterion interaction. To 
write the second term in Eq. we used Eq. Q and 
set the conditions of single slab, c = and D = 0. The 
counterion interaction with the physical charges of the 
surface is attractive and long-ranged, whereas the image- 
counterion interaction is repulsive and short-ranged. For 
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FIG. 3: (color online) The density of the counterions along 
the z direction for (a) A = 0.95 and (b) A = 0.1. In both 
plots, the solid black lines correspond to l/a = 0.3, the red 
dashed lines correspond to l/a — 0.4, and the blue dashed- 
dotted lines correspond to l/a — 0.5. 



the numerical calculations and plots, we set a = e/l 2 
and b = 2a. We note that the counterions cannot enter 
a region of thickness t on the slab. This thickness should 
be in the range of the molecular size, and we consider a 
representative value of t = 0.02a. In Fig. [2] the behav- 
ior of energy, E, is shown as a function of z for different 
values of A and I. In this figure, we set x = y = 0. As 
it can be seen, the energy has a minimum in the z di- 
rection, corresponding to a higher counterion density. In 
Fig. [3j we present the profile of the counterion density in 
the z direction for different values of A and I. The coun- 
terion interaction with the physical charges of the slab 
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FIG. 4: (color online) z cq as a function of l/a for different 
values of A. 



is attractive, whereas the image-counterion interaction is 
repulsive. One can see that by decreasing I (approaching 
to the point charge limit), the value of z corresponding 
to the maximum of the counterion density shifts toward 
zero. Decreasing of A also pushes the counterions to the 
charged surface of the slab. At very small values of I, 
independent of A, the counterions completely sit on the 
slab surface. This behavior is summarized in Fig. [4j 
where z eq is shown as a function of the charge discrete- 
ness parameter, l/a. z eq is the position of the counterion 
in the z direction, where the energy of Eq. ^ is mini- 
mum. 

It is worth saying a few words about the lateral position 
of the counterion close to the slab. When the charge 
of the slab is discretized, the counterion prefers to be 
in front of the charged regions on the slab. In Fig. [5] 
the energy of the counterion is shown in terms of x/a 
for z = z eq , A = 0.95, and different values of l/a. As 
it can be seen, the minimum of the energy corresponds 
to x/a = 0. At small values of l/a, the energy has a 
strong dependence on x/a and the depth of the energy 
well in Fig. [5]relative to ksT is quite large. The dielectric 
inhomogeneity increases the depth of the energy well and 
the tendency of the counterions to stand in front of the 
charged regions of the slab. 



Two dielectric slabs with discrete surface charge density 

Similarly, we can calculate the electrostatic energy of 
a counterion of charge q located between two dielectric 
slabs. To that end, we use Eq. [6j and find that the energy 
is minimum when the two slabs are completely in phase, 
meaning that <j> x = 4> y = (see Eq. |6]). The electrostatic 
energy of a counterion is found as 

E = q(f>2slab 
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FIG. 5: (Color online) Electrostatic energy per counterion as 
a function of x/a. The energy is calculated for z = z eg for 
each set of l/a. Plots (a) and (b) correspond to A = 0.95 and 
A = 0.1, respectively. 




where the first term comes from the electrostatic inter- 
actions between the counterion and the physical charges 
on the two slabs, and the second term takes into account 
the image-counterion interactions. In this equation, it is 
assumed that the two slabs are of the same thickness b 
and the interaction energy is found using Eq. [4] and set- 
ting c = b. By defining d = 2(D — t) and using Eq. ([7|, 
one can find the counterion density profile between the 
two slabs. It is worth mentioning that it is possible to 
suggest several values for t in the range of molecular size, 
and we consider a representative value of t = 0.02a. Fur- 
thermore, in all the numerical calculations and the plots, 
we set a — e/l 2 and b = 2a. In Fig. [6j the counterion 
density, psc, is plotted as a function of z for different 
values of A and d. Interestingly, as it can be seen in 
this figure, when the two slabs are sufficiently far from 
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FIG. 7: (Color online) The density of the counterions along 
the z direction for d/a = 2 and (a) A = 0.95 and (b) A = 0.1. 
In both plots, the solid blue lines correspond to l/a — 0.3, 
the red dashed lines correspond to l/a = 0.4, and the black 
dashed-dotted lines correspond to l/a = 0.8. 



FIG. 6: (Color online) The density of the counterions along 
the z direction for l/a = 0.2 and (a) A = 0.95 and (b) A = 
0.1. In both plots, the solid blue lines correspond to d/a — 0.1, 
the red dashed lines correspond to d/a — 0.2, the black dotted 
lines correspond to d/a = 0.4, and the green dashed-dotted 
lines correspond to d/a = 1. 

each other (for example when d > 0.15a), the counteri- 
ons prefer to be distributed in the vicinity of the slabs 
surfaces. This behavior means that the electrostatic in- 
teraction between a counterion and physical charges on 
the slabs is dominant. However, when the separation 
between the two slabs is small, the interaction between 
the counterion and the image charges dominates and the 
counterion moves to the middle of intersurface space be- 
tween the slabs. 

It is also interesting to study the effect of the discrete- 
ness of charge on the slabs. In Fig. [7] the dependence 




FIG. 8: (Color online) z eq as a function of d/a for A = 0.95 
and different values of I. 
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of psc on the parameter of the charge discreteness, I J a, 
is shown for A = 0.95 and d/a = 2. As it can be seen, 
the density of the counterions in the vicinity of the slabs 
increases with decreasing I /a. Furthermore, when the 
discreteness of charge on the slabs is washed out, the 
counterions prefer to be located in the middle of the space 
between the two slabs. To study the effects of discrete- 
ness of charge on the slabs and the distance between the 
two slabs on the preferred positions of the counterions, 
we define z eq as the preferred position of the counterion 
in the z direction, which corresponds to the position that 
p is maximum. In Fig. [8] the behavior of z eq /d is shown 
for A = 0.95 and different values of I /a. The plots show 
that for sufficiently large values of d, z eq increases as l/a 
decreases. In very small values of I, the counterions dis- 
tribute near the surfaces of the slabs at all separations. 



Interaction between the two slabs 



width of the energy wells, when the temperature is low, 
one can imagine that the counterions cannot leave their 
equilibrium positions and their fluctuations around their 
equilibrium positions is negligible. In this condition, the 
problem is equivalent to two charged slabs, with new ef- 
fective charge density, with no counterions between them. 
Due to the location of the divalent counterions on each 
slab, the effective pattern of charges of each slab follows 
a periodic positive/negative scheme on a 2D lattice of 
spacing a, as shown in the inset of Fig. [9] Therefore, the 
charge density of each slab can be written as 

Psiabi = eS(z + D) x 

+00 

(-l) m + n S(x - ma)6(y - na) 

m.n— — oo 

Psiab2 = eS(z - D) X 

+ 00 

(-l) m + n S(x -ma- <t> x )S(y - na - <j> v ), 



In this section, we study the interaction between the 
two slabs in two conditions: The limit of I — > and low 
temperatures, and the SC regime at finite temperatures. 
In the first approach despite the second one, interaction 
between the counterions is considered. In the SC regime, 
effect of the interaction between the counterions can be 
ignored. In this regime, the counterions freeze in a 2D 
lattice parallel to the surfaces, due to their strong inter- 
actions with each other relative to the thermal fluctu- 
ations 16J. The two approaches below are only appli- 
cable for divalent counterions, when counterions charge 
is twice the charge of each square-shaped region on the 
slabs (q = 2al 2 ). 



Interaction between the two slabs in the limit of I 
low temperatures 



0, and 



As one can see in Figs. [7]and[8j for very small values of 
I the counterions prefer to locate very close to the slabs. 
Therefore when the temperature is low, the entropy is 
negligible and half of the counterions are distributed in 
the vicinity of each slab. Since in this limit the electro- 
static interaction is dominant, the counterions follow the 
pattern of the slabs surface charge and make a 2D lattice 
to minimize the free energy (or the dominant electro- 
static energy). In this case the valence of the counterions 
is important. For example for mono-valent counterions, 
the 2D lattice of the counterions is exactly the same as 
the lattice of the charged regions on the slabs, whereas 
the lattice constant for the case of divalent counterions 
is y/2a. In this section, we consider divalent counteri- 
ons and find the interaction between the two slabs and 
investigate the effect of the dielectric inhomogeneity. Re- 
garding the large energy barrier against the counterions 
movement out of their equilibrium positions and small 



(10) 

where <j) x and <f) y are the relative displacement of the 
two slabs in the x and y directions. Using this charge 
distribution and the Green's function of the system, Eq. 
([2]), the electric potential can be calculated. Substituting 
this potential into the relation of energy, U = \ J p(j)dV, 
the total energy per lattice unit is found as 
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<f> x and 
the above energy 

= —1 and the total 



energy per lattice unit becomes 
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Using the relation P = — ^2 gjj" the pressure acting 



on each slab can be found as 

=2 ' J 



P 



E 

m.n— 1 



W(kb)} 2 e 2kD 



= 2fcD 



W(kb)Y 



(13) 



We see that the pressure is always negative meaning that 
the interaction is attractive. In Fig. |9j the dependence 
of pressure, P, on the separation between the two slabs, 
d, is shown for different values of A. One can see that 
by increasing the value of A, the attraction between the 
slabs becomes stronger. 
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FIG. 9: (Color online) The pressure, P, acting on each slab 
as a function of the separation between them, d. Inset: The 
pattern of the charges of each slab for divalent counterions in 
the limit of I — > and at low temperatures (see the text). 



Interaction between the two slabs in the SC regime at finite 
temperatures 

At room temperature, the total free energy of the sys- 
tem consists of the interaction energy of the two slabs 
with each other and the counterions free energy. The 
two slabs interaction energy can be obtained using the 
Green function and the slabs charge distribution as 
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In the SC limit, the counterions are frozen in a 2D lattice 
due to their strong interactions 16J. In this regime it 
can be assumed that the counterions are positioned on a 
surface parallel to the slabs at the coordinates of x = na 
and y = ma and can only move in the z direction. With 
these considerations, the system total free energy can be 
written as 



FIG. 10: (Color online) The force per unit area on a slab as a 
function of d, the distance between the two slabs, for different 
values of I. Figures (a) and (b) correspond to A = 0.95 and 
A = 0, respectively. 
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where E denotes the electrostatic energy of a counterion 
that is given by Eq. The pressure acting on a slab 
can be obtained using the relation P 



l 

2a- 



OF 
3D- 



In Fig. 10 [&) , the dependence of P, the force per unit 
area, is shown as a function of distance between the two 
slabs, d, for A — 0.95 and different values of I /a. Three 
interactions determine the total force on a slab: The 



counterions interaction with the physical charges on the 
slabs, image-countcrion interaction and the interaction 
between the physical charges of the two slabs. As one can 
see easily, the first force is always attractive, whereas the 
second and the third ones are repulsive. When the slabs 
are far from each other, the first force is always dominant 
and the total force is attractive. In closer separations, the 
effects arisen from the dielectric inhomogeneity become 
important and the force between the slabs becomes repul- 
sive. Hence, for larger values of I, an equilibrium separa- 
tion between the slabs exists. This is in agreement with 
the result of Ref. [14] for uniform distribution of charge 
on the slabs. But for smaller values of /, when the slabs 
are in close separation, the electrostatic attraction dom- 
inates again and the slabs completely attract each other 
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after passing a repulsion barrier. It is interesting to note 
that the latter behavior is completely different from the 
behavior of the slabs with uniform surface charge which 
shows the important role of the charge discreteness. 

The behavior of P is shown in Fig. [To|b) as a function 
of d for A = and different values of I /a. As one can 
see, when I is large and the charge distribution on the 
slabs is approximately uniform, the force on the slabs in 
close separations is repulsive which is the result of the 
counterions entropy |17j . The interaction force becomes 
attractive for large separations and there is an equilib- 
rium distance for the two slabs. Furthermore, for small 
values of I, the electrostatic attraction becomes stronger 
and overcomes the repulsion. As it can be seen in this 
figure, at small values of I the two slabs attract each other 
in all separations. This is in agreement with the result 
of Ref. [115] that the charge discreteness strengthens the 
attraction between the slabs. 

In Figs. |ll[ a)— (c), attraction and repulsion regions in 
A-d plane are shown for the system. One can see that 
by increasing the charge discreteness, the structure of 
the phase space changes completely. By decreasing I, the 
repulsion region in small separations disappears and a 
repulsion region appears in the intermediate separations. 
At each value of I, with increasing A the repulsion re- 
gion grows. The interesting point in these plots is that 
for specific values of I and A, there are two equilibrium 
separations between the two slabs. 

In Fig. [12] the behavior of the free energy of the system 
is shown as a function of d/a for different values of l/a 
and A = 0.8. As one can see, for the values of l/a > 0.25, 
the free energy has two minima showing two equilibrium 
separations between the two slabs. It is interesting to 
note that these two minima have the same free energy 
for l/a ~ 0.195. For l/a < 0.195, the free energy has 
a global minimum when the distance of the two slabs is 
about d/a < 0.2. When l/a > 0.195, the global minimum 
of the free energy tends to the larger distance of the two 
slabs, say d/a > 0.5. As it can be seen in the figure, for 
the larger values of l/a (e.g. l/a > 0.26), the free energy 
has only one minimum that is located in d/a > 0.5. 

DISCUSSION AND CONCLUDING REMARKS 

Strongly-charged biological macromolecules such as 
DNA and filamentous actin are in the limit of strong cou- 
pling in biological conditions. For these macromolecules 
in solution, it is known that multivalent counterions are 
mostly distributed in their vicinity. When the small ions 
are very close to these macromolecules or when two of 
them are very close to each other, the curvature of their 
surface can be ignored in calculation of the electrostatic 
interactions. Therefore, our study of the interaction be- 
tween two flat dielectric slabs could be useful in enlight- 
ening the phenomenon of interaction between biological 



charged macromolecules in the presence of multivalent 
counterions. As examples, like-charge attraction between 
actin filaments [25], the interaction between DNA and 
charged proteins |2. and a system of proteins confined in 
a polyelectrolyte brush J35] can be considered. In most 
of computer simulation studies of such systems dielectric 
inhomogeneity and real charge distribution on the surface 
of the macromolecules have not been taken into account. 
This is because of intrinsic difficulty of calculation of the 
electrostatic interactions in the presence of dielectric in- 
homogeneities and massive calculations needed for mod- 
eling the real charge distribution on the macromolecules. 
In this subject only simple systems like single or a cou- 
ple of dielectric slabs with uniform charge distribution on 
their surfaces have been studied. 

By calculation of the Green function for two dielec- 
tric slabs (the geometry shown in Fig. [I]), we calcu- 
lated the electrostatic interaction between them taking 
into consideration both the dielectric inhomogeneity and 
the charge discreteness. In the strong coupling limit, it 
is found that the amount of dielectric constant difference 
between the slabs and the environment and the discrete- 
ness of charge on the slabs have opposing effects on the 
equilibrium distribution of the counterions between the 
slabs. Increasing the amount of dielectric constant dif- 
ference increases the tendency of the counterions toward 
the middle of the space between the slabs. Discreteness 
of charge on the slabs however, pushes the counterions 
to the surface of the slabs. At low temperatures, the in- 
teraction between the slabs is attractive, the equilibrium 
separation between the slabs vanishes and the strength 
of their attraction increases with increasing the amount 
of the dielectric constant difference. At room tempera- 
ture, the slabs may attract each other and come together, 
stay in an equilibrium separation or have two equilibrium 
separations with a barrier in between depending on the 
system parameters. 

Our results showed that for a system of dielectrics 
of a simple geometry, namely two dielectric slabs, tak- 
ing into consideration the dielectric inhomogeneity and 
leaving the reality of discreteness of charge on the di- 
electrics is not a valid approach when they are close to 
each other. In the extreme of point charges on the di- 
electric slabs, the effect of charge discreteness completely 
dominates over the effect of dielectric inhomogeneity. For 
example, distribution of counterions in the middle of the 
two slabs because of the dielectric inhomogeneity (as de- 
scribed in Rcf. [14]) is the effect of the assumption of 
uniform charge distribution. Our results show that by 
taking into account the reality of charge discreteness on 
the slabs, all the counterions distribute in the vicinity of 
the slabs surfaces and the middle of the slabs correspond 
to the minimum of the counterions density. 

When the surface charge density is uniform, the pres- 
sure between the plates is repulsive at small separations 
and the plates stand in an equilibrium distance from each 
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FIG. 11: (Color online) The phase space of the system for three representative cases of (a) / = 0.01a, (b) / = 0.15a, and (c) 
I — 0.7a. The plots are delineating the different regimes in the parameter space of two charged slabs, where d and A denote the 
distance between the two slabs and dielectric inhomogeneity, respectively. The solid line corresponds to P = and the dashed 
line corresponds to dP/dA — 0. 
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FIG. 12: (Color online) Free energy of the system for A = 0.8 
and several values of I. At l/a ~ 0.195, two minima of the 
free energy have the same value (the horizontal dashed line). 



other. This repulsion is due to the counterions entropy 
and the dielectric inhomogeneity. As the charge discrete- 
ness increases, the repulsion at close separations disap- 
pears and the plates completely attract each other after 
passing a repulsion barrier in the intermediate separa- 
tions. Also, dielectric inhomogeneity widens the repul- 
sion region between the slabs and increases the repulsion 
strength. For highly discretized surface charge, the inter- 
action between the slabs was shown that is attractive at 
small inter-slab separations. In this case, the dielectric 
inhomogeneity may increase or decrease the attractive 
pressure, depending on the value of the parameters. 

Another point to note is the large depth of the en- 
ergy wells in Figs. [2] and [5] relative to the thermal en- 
ergy, ksT. It shows that the charge discreteness and the 



dielectric inhomogeneity strongly affect the equilibrium 
distribution of the counterions and hence the interaction 
between the slabs. Regarding the competing effects of 
charge discreteness and the dielectric inhomogeneity, our 
results showed that considering only one of these reali- 
ties and leaving out the other one may crucially bias the 
results of the theory. 

In this paper, distribution of the counterions and in- 
teraction between the slabs are obtained from the SC 
theory at the leading order of the expansion. It has been 
shown that this approximation is correct for discretely 
charged surfaces in coupling parameters larger than 20 
|19j . Many biological systems lie in this range. However, 
when the separation between the two charged plates gets 
very large, the first order term becomes smaller than the 
second term and this approximation fails. This critical 
distance is found for uniformly charged plates with and 
without dielectric inhomogeneity. Computer simulations 
have also been used to find the validity range of the SC 
theory in these systems [THEI]. As a result, computer 
simulations considering both the dielectric inhomogene- 
ity and the charge discreteness effects seem essential. 

The SC theory is obtained from an expansion of the 
system parameters in inverse powers of the coupling pa- 
rameter 3. Theoretically, SC formalism is exact in the 
limit of 3 — > oo, however computer simulations have 
shown that the results become accurate for values of 
5 > 100 [IB]. In biological conditions, 3 = 100 is equiv- 
alent to a ~ lnm in our model or surface charge density 
as — 1 which is of the order of the typical charge den- 
sity of charged biological membranes and DNA [17]. For 
discretely charged surfaces, SC theory has been shown 
to be exact even in lower values of 5 (down to 3 ~ 20) 
due to the weakened correlations between the counteri- 
ons [TH] • Also, it has been shown that the presence of the 
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dielectric inhomogeneity widens the validity range of the 
SC theory [14]. 

The validity range of the SC theory in the case of a 
single charged plate is z < a± in which z is the sep- 
aration between counterions and the surface and a± is 
the mutual separation between the counterions parallel 
to the surface. This range for the system of two parallel 
charged surfaces is d < a± in which d is the separation 
between the surfaces [IB] • In our model, the lattice spac- 
ing of the regular charge distribution on the surfaces, a, 
plays the role of a± of the SC theory. Considering that 
the most of our results have been obtained using param- 
eters in the range of d < a, the validity of the SC theory 
is preserved. 

One should note that a sharp boundary is assumed be- 
tween inside and outside of the slabs in our model. In 
reality, it is known that such a sharp boundary is not the 
case and an improvement of our model could be the con- 
sideration of a smooth profile of the dielectric constant 
in the slabs boundaries. Some additional points in the 
way of approving the theory are as follows. The grand 
canonical ensemble is more proper for study of the real 
systems such as the system of dielectric slabs with coun- 
terions studied here but the SC theory is developed in 
the canonical ensemble. Developing a theory in which 
the number of the counterions between the slabs is not 
constant seems as a valuable step. The finite distance of 
the fixed charges from the slab surface [TU [T!5] , random 
distribution of the charged domains [27] , displacement of 
the charged domains by changing the system parameters 
in some cases [35J and possibility of their hydration [21] 
should be considered in development of the theory. Also, 
experimental studies on the effects of different factors on 
the distribution of the counterions and interaction be- 
tween charged dielectrics can also be insightful. For ex- 
ample, the effect of the dielectric inhomogeneity on the 
interaction of two charged dielectrics seems to be possi- 
ble by changing the solvent [3D] and is not studied yet to 
the authors knowledge. 



ACKNOWLEDGMENT 

We are grateful to F. Julicher, S. N. Rasuli, and M. F. 
Miri for interesting discussions and comments. 



* Electronic address: farshid@iasbs.ac.ir 
' Electronic address: fazli@iasbs.ac.ir 
[1] O. Halskau et al. , Proc. Natl. Acad. Sci. U.S.A. 105, 
8625 (2008); D. R. Ripoll, J. A. Vila, and H. A. Scheraga, 
Proc. Natl. Acad. Sci. U.S.A. 102, 7559 (2005); F. Avbelj 
and R. L. Baldwin, Proc. Natl. Acad. Sci. U.S.A. 100, 
5742 (2003). 



[2] V. Dahirel et al., Phys. Rev. Lett. 102, 228101 (2009); D. 
Jiao, P. A. Golubkov, T. A. Darden, and P. Ren, Proc. 
Natl. Acad. Sci. U.S.A. 105, 6290 (2008); P. Gong and R. 
Levicky, Proc. Natl. Acad. Sci. U.S.A. 105, 5301 (2008); 
M. Borg et al., Proc. Natl. Acad. Sci. U.S.A. 104, 9650 
(2007). 

[3] M. K. Krotova et al, Phys. Rev. Lett. 105, 128302 
(2010); S. Lee, T. T. Le, and T. T. Nguyen, Phys. Rev. 
Lett. 105, 248101 (2010); X. Qiu et al., Phys. Rev. Lett. 
101, 228101 (2008); V. A. Belyi and M. Muthukumar, 
Proc. Natl. Acad. Sci. U.S.A. 103, 17174 (2006); N. V. 
Hud and K. H. Downing, Proc. Natl. Acad. Sci. U.S.A. 

98, 14925 (2001). 

[4] T. S. Lo, B. Khusid, and J. Koplik ,Phys. Rev. Lett. 

100, 128301 (2008); P. Gong, J. Genzer, and I. Szleifer, 

Phys. Rev. Lett. 98, 018302 (2007); Y. Mei et al., Phys. 

Rev. Lett. 97, 158301 (2006); P-.Y. Hsiao and E. Luijten, 

Phys. Rev. Lett. 97, 148301 (2006). 
[5] M. Tagliazucchi, M. O. de la Cruz, and I. Szleifer , Proc. 

Natl. Acad. Sci. U.S.A. 107, 15300 (2010); P. Besenius 

et al, Proc. Natl. Acad. Sci. U.S.A. 107, 17888 (2010); 

G. H. Lai et al, Phys. Rev. Lett. 101, 218303 (2008); M. 

O. Steinmetz et al, Proc. Natl. Acad. Sci. U.S.A. 104, 

7062 (2007); L. K. Sanders et al, Proc. Natl. Acad. Sci. 

U.S.A. 104, 15994 (2007); H. Liang, D. Harries, and G. 

C. L. Wong, Proc. Natl. Acad. Sci. U.S.A. 102, 11173 
(2005); D. J. Needleman et al, Proc. Natl. Acad. Sci. 
U.S.A. 101, 16099 (2004). 

[6] F. Huang et al, Phys. Rev. Lett. 102, 108302 (2009); D. 
K. Satapathy et al, Phys. Rev. Lett. 101, 136103 (2008); 
A. -P. Hynninen and A. Z. Panagiotopoulos, Phys. Rev. 
Lett. 98, 198301 (2007). 

[7] R. Das et al., Phys. Rev. Lett. 90, 188103 (2003); X. Qiu 
et al, Phys. Rev. Lett. 96, 138101 (2006). 

[8] A. A. Kornyshev and S. Leikin, Phys. Rev. Lett. 82, 
4138 (1999); X. Qiu et al, Phys. Rev. Lett. 99, 038104 
(2007); L. Dai, Y. Mu, L. Nordenskiold, and J. R. C. van 
der Maarel, Phys. Rev. Lett. 100, 118301 (2008); X. Qiu, 
V. A. Parsegian, and D. C. Rau, Proc. Natl. Acad. Sci. 
U.S.A. 107, 21482 (2010). 

[9] E. E. Meyer et al., Proc. Natl. Acad. Sci. U.S.A. 102, 
6839 (2005); M. Trulsson et al., Phys. Rev. Lett. 97, 
068302 (2006); O. Zohar, I. Leizerson, and U. Sivan, 
Phys. Rev. Lett. 96, 177802 (2006); E. Greenfield and 
U. Sivan, Phys. Rev. Lett. 102, 106101 (2009). 
[10] J. A. Freites, D. J. Tobias, G. von Heijne, and S. H. 
White, Proc. Natl. Acad. Sci. U.S.A. 102, 15059 (2005); 
V. Gonzlez-Prez, K. Stack, K. Boric, and D. Naranjo, 
Proc. Natl. Acad. Sci. U.S.A. 107, 5178 (2010); S. 
Buyukdagli, M. Manghi, and J. Palmeri, Phys. Rev. Lett. 
105, 158103 (2010). 
[11] H. I. Petrache, T. Zemb, L. Belloni, and V. A. Parsegian, 
Proc. Natl. Acad. Sci. 103, 7982 (2006); J. Zwanikken 
and R. van Roij, Phys. Rev. Lett. 99, 178301 (2007); 

D. Frydel, S. Dietrich, and M. Oettel, Phys. Rev. Lett. 

99, 118302 (2007); Y. Levin, Phys. Rev. Lett. 102, 
147803 (2009); K. Masschaele et al, Phys. Rev. Lett. 
105, 048303 (2010). 

[12] J.-Y. Wang et al, Phys. Rev. Lett. 96, 128301 (2006). 
[13] P. Sens and J.-F. Joanny, Phys. Rev. Lett. 84, 4862 
(2000). 

[14] Y. S. Jho et al, Phys. Rev. Lett. 101, 188101 (2008). 
[15] Y. S. Jho et al, Phys. Rev. E 76, 011920 (2007). 
[16] R.R. Netz, Eur. Phys. J. E 5, 557 (2001). 



13 



[17] A. Naji, S. Jungblut, A. G. Moreira, R. R. Netz, Phys. 

A 352, 131 (2005). 
[18] A. G. Moreira and R. R. Netz, Europhys. Lett. 57, 911 

(2002). 

[19] Y. S. Jho et al, Phys. Rev. E 73, 021502 (2006). 
[20] A. Naji et al, Phys. Rev. Lett. 104, 060601 (2010). 
[21] S. Buyukdagli, M. Manghi, and J. Palmeri, Phys. Rev. E 

81, 041601 (2010); M. M. Hatlo and L. Lue, Europhys. 

Lett. 89, 25002 (2010). 
[22] L. Samaj, and E. Trizac, Phys. Rev. Lett. 106, 078301 

(2011). 

[23] D. B. Lukatsky, S. A. Safran, A. W. C. Lau, and P. Pin- 
cus, Europhys. Lett. 58, 785 (2002); M. L. Henle, C. 
D. Santangelo, D. M. Patel and P. A. Pincus, Europhys. 
Lett. 66, 284 (2004). 



[24] D. B. Lukatsky and S. A. Safran, Europhys. Lett. 60, 
629 (2002). 

[25] G. C. L. Wong et al, Phys. Rev. Lett. 91, 018103 (2003). 
[26] K. Henzler et al, Phys. Rev. Lett. 100, 158301 (2008); 

K. Henzler et al, J. Am. Chem. Soc. 132, 3159 (2010); 

K. Henzler et al, Phys. Chem. Chem. Phys. 13, 17599 

(2011). 

[27] A. Naji et al, Phys. Rev. Lett. 104, 060601 (2010). 
[28] R. Brewster, P. A. Pincus, and S. A. Safran, Phys. Rev. 

Lett. 101, 128101 (2008). 
[29] D. Hollenbeck et al, Phys Rev. E 82, 031402 (2010). 
[30] K. Besteman, M. A. G. Zevenbergen, and S. G. Lemay, 

Phys. Rev. E 72, 061501 (2005). 



